How accesss to care shapes disease burden: the current impact of post-exposure prophylaxis and potential for expanded access to prevent human rabies deaths in Madagascar
Malavika Rajeev1,✉, Helene Guis4,6, Glenn Edosoa7, Chantal Hanitriniaina8, Soa Fy Andriamandimby2, Anjasoa4, Reziky4, Ravo3, Laurence Baril5, C.J.E. Metcalf1, and Katie Hampson9
1 Department of Ecology and Evolutionary Biology, Princeton University, Princeton, United States
2 Virology Unit, Institut Pasteur de Madagascar, Antananarivo, Madagascar
3 Vaccination Center, Institut Pasteur de Madagascar, Antananarivo, Madagascar
4 Epidemiology and Clinical Research Unit, Institut Pasteur de Madagascar, Antananarivo, Madagascar
5 Institut Pasteur de Cambodge
6 CIRAD, UMR ASTRE, Antananarivo, Madagascar
7 OMS
8 Mention Zoologie et Biodiversité Animale, Faculté des Sciences, Université d’Antananarivo, Antananarivo, Madagascar
9 Boyd Orr Centre for Population and Ecosystem Health, Institute of Biodiversity, Animal Health and Comparative Medicine, University of Glasgow, Glasgow, UK
✉ Correspondence: Malavika Rajeev <mrajeev@princeton.edu>
ABSTRACT
INTRODUCTION
METHODS
Spatial covariate data
To calculate travel times, we used the global friction surface for 2015 generated by the Malaria Atlas Project (https://map.ox.ac.uk/research-project/accessibility_to_cities/, [1]) and GPS points of clinics to estimate the travel time to the nearest ARMC for the country of Madagascar at an approximately 1 x 1 km scale. We calculated a weighted average of travel times by human population at the commune and district levels, using administrative shapefiles from the UN Office for the Coordination of Humanitarian Affairs (see section S1). Human population estimates were taken from the 2015 UN adjusted population projections from World Pop (www.worldpop.org, [2]) and aggregated to the commune and district levels.
Bite patient data
We used a database of bite patient forms submitted to IPM from ARMC across the country between 2014 - 2017. These were individual patient data forms that were submitted from clinics as frequently as monthly to annually. These data include details of the administrative district of the bite patient’s address and the date of reporting.
For most districts, the majority of bites were reported to the closest clinic as estimated by our travel time metric (Fig 1A, Fig S1.2). Therefore, we defined the catchment for each clinic as all districts for which the clinic was the closest ARMC for the majority of the population in that district. We excluded any districts in a catchment of a clinic which submitted less than 10 forms (excluded NN catchments, corresponding to XX districts, in grey in Fig 1A).
FIG1 A caption + fig should go here!
Figure 1. A) The network of patient presentations at the district level for the national data: circles with black outline represent the total number of patients reporting to each clinic for which we have data. Color corresponds to the clinic catchment. Each filled circle is the total number of bites reported for that district, and the lines show which ARMC those patients reported to, with the line width proportional to number of patients. Districts in catchments excluded due to lack of forms submitted by the clinic are colored in grey. B) The estimated annual bite incidence corrected for reporting for the years 2014 - 2017 (some years were excluded due to low reporting) for each district with colors corresponding to the catchment and the y-axis ordered by district level travel times to the nearest ARMC. C) The network of patient presentations at the commune level for the Moramanga data with inset of Madagascar showing the location of the district. D) The estimated annual bite incidence for each commune in the Moramanga catchment with the y-axis ordered by commune level travel times.
For clinics that did submit data, there was substantial underreporting, with months to even a whole year when forms were not submitted. We therefore estimated clinic level reporting as the proportion of days on which forms were submitted, excluding any periods for which there were no forms submitted for 15 consecutive days (Figure S2.1A and B). To estimate the average annual bites reported for each district, we further excluded data from any years for which there was less than 25% reporting at the clinic level. We found that these corrected estimates of reported bites was in agreement with the vial demand for clinics (Fig S2.1C, further details in section S2).
Previous work in the Moramanga District showed that Category I exposures (cite WHO), i.e. touching or feeding a suspected animal or human rabies case, make up approximately 20% of patients reporting to ARMCs [3]. These low risk exposures do not contribute to the mortality burden of rabies, however we could not directly exclude them as details on exposure type was missing for the majority of the patient forms. However, detailed data from Moramanga indicates that these contacts consistently present as clustered cases. Accordingly, to reflect this observation, we excluded patients reporting on any dates that had greater than 3 standard deviations above the mean number of patients reporting per day (Figure S2.2A). We validated this method using the Moramanga ARMC data for which we had details on the type of exposure, and found that setting the threshold to 3 standard deviations (SDs) excluded approximately 50% of known contacts, with only 2% of Category II and III patients excluded (Figure S2.2B). For the national data for which a subset of patient forms were explicitly noted to be contacts, we found that our exclusion criteria of 3 SDs identified 68.28% of known contacts. We further excluded known contacts which were not identified based on this method, resulting in the exclusion of approximately YY.Y % of patient records from the national data.
After excluding patients likely to reflect low risk contacts as described, we estimated average bite incidence as the mean of the yearly bites reported/the proportion of bite patient data which were submitted. Our final dataset consisted of estimates of average bite incidence for 81 of 114 districts from 20 catchments (Figure 1A and 1B).
The 31 months of data available from the Moramanga district resolved to the commune level (the administrative level below the district, Figure 1C) provide a useful comparison to the national scale data. These data were completely reported with details on the type of exposure documented. We calculated average annual bite incidence to the commune level by excluding Category I exposures and taking the monthly average of reported bites and multiplying by 12 (Figure 1D).
Model of reported bites as a function of travel time to closest ARMC
We modeled incidence of bites as a function of travel time (\(T\)) as follows: \[ \mu = exp(\beta_{t}T \ + \ \beta_0)P\] where \(\mu\) is the expected number of bite patients presenting at the ARMC as a function of travel time to the ARMC, and the human population size, \(P\) (an offset which scales the predicted number of bites), for a given source location (district or commune). We fit this model to both the national data resolved to the district level and the Moramanga data resolved to the commune level in a Bayesian framework via MCMC using the R package ‘rjags’. In order to more directly compare estimates between these two datasets, we also modeled the national data with a commune-level travel time covariate (\(T_{j}\)) so that: \[\mu = \sum \limits_{j=1}^jexp(\beta_{t}T_j\ + \ \beta_0)P_j\]
As travel times are correlated with human population size (Figure S3.1A), we also compared how well bites were predicted by human population size alone, and in combination with travel time. For the models with population size, we used a model framework with either population size alone (\(\mu = exp(\beta_pP + \beta_0)\)) or population size and travel times (\(\mu = exp(\beta_{t}T + \beta_pP + \beta_0)\)) as predictors of bites for the Moramanga data, the national data with covariates at the district level, and the national data with covariates at the commune level.
For the models fitted to the national data, we also attempted to capture variation between clinics by modeling a random catchment effect (\(B_0_k = norm(\mu_k, \sigma_k\)) where \(\mu_j\) is the mean and \(\sigma_j\) is standard deviation for the intercept for each catchment \(k\)). To test our models predictions further, we used parameter estimates from each model to generate predictions for data that were not used to fit the model (i.e. we used the estimates from models fitted to the national data to predict the Moramanga data, and used estimates from models fitted to the Moramanga data to predict the national data).
Estimating burden of rabies deaths
We estimate deaths due to rabies as a function of the number of bites predicted by our model and estimates of endemic rabies exposure incidence using an adapted decision tree framework (Fig 2). We assume annual rabies incidence of 1% in the dog population in an endemic setting with no mass dog vaccination (cite GAVI). To generate \(E_i\), the human exposure incidence in a given administrative unit (where \(i\) is the index of the commune or district) based on the model, we multiply this baseline incidence by the estimated dog population and the exposure rate per rabid dog (\(p_{exp}\)). We estimate the dog population given a human:dog ratio (HDR) of 7 to generate our maximum expected incidence and an HDR of 35 for our minimum expected incidence, giving a range of 11.1 (\(E_{min}\)) - 55.7 (\(E_{max}\)) exposures per 100,000 persons. As there is little data on dog populations in Madagascar, this range encompasses both a wide range of observed HDRs and human exposure incidence across Africa (cite). \(E_i\) is then drawn from a uniform distribution between the minimum and maximum expected number of human exposures for each location.
Figure 2.Decision tree for burden estimation. For a given administrative unit \(i\), human deaths due to rabies (\(D_i\)) are calculated from model predicted bites (\(B_i\)). To get \(R_i\), the number of reported bites that were rabies exposures, we multiply \(B_i\) by \(p_{rabid}\), the proportion of reported bites that are rabies exposures. \(p_{rabid}\) is constrained such that reported rabies exposures cannot exceed total rabies exposures (\(E_i\)) and reporting cannot exceed a maximum (\(p_{max}\)). \(E_i\) is generated from assumptions about rabies exposure incidence and the underlying dog population extrapolated from human population sizes (see Table 1, Methods text). \(R_i\) is subtracted from \(E_i\) to get the number of unreported bites (\(U_i\)) and then multiplied by the probability of death given a rabies exposure (\(p_{infect}\)) to get the number of deaths due to rabies (\(D_i\)). Similarly, deaths averted by PEP are estimated by multiplying \(R_i\) by \(p_{infect}\), i.e. those who would have died given exposure, but instead received PEP. All parameter values and sources are listed in Table 1.
\(B_i\), the number of total reported bites, is drawn from a poisson distribution with the mean predicted number of bites predicted from the models of incidence as a function of travel times. \(p_{rabid}\) is the proportion of the reported bites which we assume are rabid exposures. We draw \(p_{rabid}\) from a uniform distribution ranging from 0.2 - 0.6 as reported in Rajeev et al. 2018 based on investigations of bite patients reporting in the Moramanga District. So that rabid reported bites cannot exceed the total expected number of human rabies exposures or a maximum reporting (as we know that even with minimal travel times, people may not report for PEP for other reasons), we constrain \(p_{rabid}\) as follows: \[ p_{rabid}=\begin{cases} x, & \text{if}\ \frac{E_i\rho_{max}}{B_i} > x \\ \frac{E_i\rho_{max}}{B_i}, &\text{otherwise} \end{cases}\] where \(\rho_{max}\) is the maximum reporting. We estimate \(\rho_{max}\) from the Moramanga ARMC data for the commune of Moramanga Ville, the closest commune to the ARMC (average of 3.12 minutes travel time to the clinic), and where contact tracing data estimated that approximately 2% of rabies exposure go unreported (Rajeev, unpublished data). We assume that all rabies exposed patients who report to an ARMC receive and complete PEP, and PEP is completely effective at preventing death due to rabies. We subtract \(B_i\) from \(E_i\) to get the number of unreported exposures (\(U_i\)). Deaths due to rabies in each administrative unit (\(D_i\)) are then estimated as: \[ D_i = binom(U_i,\ p_{death}) \] Deaths averted are calculated similarly as the reported rabies exposures \(R_i\) that would have died in the absence of PEP given \(p_{death}\). Table 1 lists all parameter values and sources used in the decision tree.
| Parameter | Value | Description | Source |
|---|---|---|---|
| \(p_{exp}\) | 0.39 | Probability of human exposure given a dog rabies case | [4] |
| \(E_{min} − E_{max}\) | 11.1 - 55.7 per 100,000 persons | Range of human exposure incidence calculated from an assumed dog rabies incidence of 1%, human:dog ratios between 5 - 25, and \(p_{exp}\). | [4–6] |
| \(p_{rabid}\) | 0.2 - 0.6 | Proportion of reported bites that are rabies exposures | [3] |
| \(\rho_{max}\) | 0.98 | The maximum reporting possible for any location Data from Moramanga Ville, Moramanga District | [3] |
| \(p_{death}\) | 0.16 | The probability of death given a rabies exposure | [7] |
Strategically expanding PEP access
We developed a framework to rank clinics by how much they could improve access for people currently underserved (< 3 hours away from any clinic in the country) and to estimate the incremental reduction in burden as PEP is expanded to these clinics. We use the approximate locations of all existing public CSB II (centre de sante niveau II, the level below district hospitals, n = 1679) as locations of potential ARMC.
Starting with the current locations, we added one clinic at a time, calculating the proportion of people living < 3 hours away from any clinic for the country. We added the clinic which minimized this metric weighted by the reduction in travel times, and then repeated the process iteratively, ranking clinics and adding the top clinic sequentially. In this iterative process, we also exclude any clinics that once added only change travel times above the 3 hour threshold for less than 0.01% of the population (~2300 persons). We also estimate travel times given that all CSB II in the country provision PEP to get an estimate of the maximum impact. Finally, we estimate the incremental reduction in burden as each clinic is added using the burden framework described in the previous section.
We also estimate how expansion to additional clinics may impact vaccine demand. Specifically, we use our model-predicted estimates of annual bites and then aggregate these to the clinic catchment (where the catchment clinic is the one which is the closest for the majority of people living in the district). For the district model, we disaggregate the reported bites to the commune level (assuming uniform incidence within the district) so as to better assign bites to a given catchment (as multiple clinics now serve a single district, see Figure S5.1?).
For each catchment, we simulate throughput by randomly assigning patient presentation dates, and then assume perfect compliance (i.e. all patients reporting subsequently on day 3 and day 7 for their second and third doses) to generate dates of subsequent doses. We use these dates to estimate vial usage given routine vial sharing practices in Madagascar (i.e. each vial is split between two patients, and vials are discarded at the end of each day) and the new abridged intradermal regimen recommended by WHO (2 injection sites on day 0, 3, and 7, [8]). We assume that vial sharing is functioning perfectly, with a background wastage of 0.1 mL per vial, but not accounting for wastage due to other reasons (i.e. wastage due to injection error). We take the mean of 100 simulations of throughput for each catchment and set of clinics.
Sensitivity analyses for burden
We estimate burden across a range of parameter values to test the effects of our model assumptions on estimates of rabies burden. Specifically, we fix rabies incidence at the minimum and maximum of our estimated range, \(p_{rabid}\) at 0.2 and 0.6, and \(\rho_{max}\) at 0.85 and 0.99 to obtain maximum and minimum estimates of burden. We also examine the impact of systematic variation in rabies incidence with human population (a potential proxy for changes in the dog population, see section S5) by looking at how estimates of burden change if rabies incidence were to scale positively or negatively with human population size.
RESULTS
Travel times as a predictor of bites
All of the candidate models which included travel times produced reasonable predictions when compared to the observed data (Figure SY). Population size alone was the poorest fit to the data as estimated by DIC, and all models with population size as a covariate did not generate realistic predictions to the data or when used to predict out of fit. For the national data, the models with travel times alone as a covariate and a random catchment generated predictions closest to the observed data (Figure 3B, Figure SX). In addition, the estimates from the model fitted to the Moramanga data fall into the prediction envelope for the random effect models (Figure 3A). As these models likely capture the range of possible relationships for any given catchment in our data set, we use estimates from these models to generate subsequent predictions of bite incidence and burden.
Figure 3. Travel times as a predictor of reported bite incidence per 100,000 persons A) the estimated relationship between travel times in hours (x-axis) and mean annual bite incidence per 100,000 persons (y-axis). The lines are the estimates using the mean of the posterior distributions of the parameters for three candidate models: model with travel times at the 1) commune and 2) district level fitted to the national data with a random intercept by catchment, and 3) travel times at the commune level fitted to the Moramanga data with a fixed intercept. The envelopes shows the 95% CI given the random intercept for the commune and district models. The Moramanga data was fitted with a fixed intercept, but falls in the prediction envelope for the commune model with a random intercept. B) Log(observed) vs. log(predicted) mean annual reported bites for the three different models (added a small positive offset of 0.1 to account for log of zeros).
Estimation of baseline burden
Overall, we estimate between 306 - 373 average annual deaths in Madagascar due to rabies depending on the model used (Table 2). We also estimate that under the current system of 31 ARMCs in Madagascar, approximately 483 - 548 deaths due to rabies are prevented through provisioning of PEP each year. In general, the incidence of rabies deaths increases with travel times to clinics, and there is significant variation sub-nationally at the district level, but also within districts (Figure 4).
Figure 4. Spatial variation in predicted incidence of human rabies deaths per 100, 000 persons. A) for each district (y-axis) in Madagascar, where the diamonds show the predicted incidence for the district model and the squares show the incidence for the commune model for all communes in a district. Points are colored by and districts are ordered by travel times. The dashed lines show the average national incidence of human rabies deaths for the commune (red) and district (blue) models. Incidence mapped to the B) commune level and C) district level from the respective models; grey x’s show locations of current CTAR.
Table 2. Model predictions of mean annual reported bite incidence, total deaths, and total deaths averted at the national level for the two models (commune level and district level models with a random catchment effect); 95% CIs in parentheses.
Incremental reduction of burden with expanded PEP
We use our travel time models to look at the benefits of expanding PEP incrementally. We find that expanding access to PEP reduces rabies deaths, but benefits accrue less rapidly as ARMC are added, saturating after approximately 200 clinics. This effect is more pronounced for the commune level model than for the district model. Vaccine demand also increases, and more vials are needed per averted death as PEP access is expanded, with vials needed outpacing the increase in deaths averted. Report actual # of vials needed as PEP is expanded. Even with maximal PEP access, our model still predicts > 100 rabies deaths per annum, as average reporting stays at x% (see section S6 to see how things shift!)
Figure 5. A) Decrease in deaths due to rabies and B) increase in vials needed per death averted as additional ARMC are added at the national level based on the two models of reported bites. Lines are the mean of 1000 simulations with envelopes representing 95% CIs. The “…” on the axes indicates a break in the axes, and the last point is the burden and vial predictions when all CSBS in the country provision PEP (max access).
Sensitivity analyses
There is significant uncertainty in the predicted numbers of rabies deaths, but the broad patterns of higher risk with increased travel times are robust for both our estimates of burden and the impacts of expanding access to PEP, even given systematic bias in our assumptions of incidence (see Supplementary Materials, section S5). Our estimate of the absolute number of deaths is most sensitive to our assumptions of rabies exposure incidence in the population.
DISCUSSION
Strengths and Limitations
Broader Context
Conclusion
SUPPLEMENTARY MATERIALS AND METHODS
S1. Estimating travel times to the nearest ARMC
jump back to main methods section
To estimate travel times to the nearest clinic at the scale of administrative units, we generated two raster inputs: 1) the travel time estimates generated using the friction surface from the Malaria Atlas Project [1] at an ~1 km^2 scale and the locations of existing ARMC and 2) the population estimates which we resampled from the World Pop (worldpop.org) raster of population originally at an ~100m^2 resolution. We extracted the mean of travel times (Fig S1.1A) weighted by the population in that grid cell (Fig S1.1B) for each district or commune. Figure S1.1 shows the rasters we generated as well as the resulting extracted estimates of travel times at the district (Fig S1.1D) and commune (Fig S1.1D) level. Raw inputs are availble through World Pop and MAP[9].
Figure S1.1. A) Travel time estimates at an approximately 1x1 km scale B) Population per 1x1 km grid cell. Mean travel times weighted by population size extracted for each C) District and D) Commune.
We assigned clinic catchments by determining which clinic was closest in terms of travel times for the majority of the population within the administrative unit. Figure S1.2A shows that for most districts and communes a single clinic served the majority of the population by this metric. In addition, our data show that the majority of bite patients do report to the closest clinic as assigned by our travel time metric (Fig S1.2B).
Figure S1.1. A) Distribution of the proportion of the population in a given administrative unit (district or commune) served by the catchment assigned. B) The proportion of bites reported to each clinic which originated from a district within the assigned catchment.
S2. Correcting data for underreporting and excluding Category I exposures
Estimating under-submission of forms
To correct for under-submission of forms, we attempted to identify what proportion of data were missing each year by excluding any periods of time for which there were 15 consecutive days where zero patient records were submitted (Figure S2.1A). We validated this cut-off of consecutive days by comparing estimated vial demand given the total reported bites (corrected for under-submission) for the years 2014-2016 to actual vials provisioned to clinics over this same period. We did not have the date or month when vials were provisioned, so we compared the totals across this three year period.
Vial demand was simulated under simplified assumptions of PEP administration and adherence[3], assuming that patients reported randomly across the year and then completed doses 3, 7, and 28 days subsequent from the date reported. For this analysis, we look at the range of vial estimates assuming that all patients complete 3 vs. 4 doses total. During this period, the Thai Red Cross Intradermal regimen was used across Madagascar, with 0.2 mL administered per patient at each visit. Vials can be shared within a day between two patients, resulting in 0.1 background wastage per vial shared + any additional wasteage from unused doses discarded at the end of the day. While these simplifying assumptions do not capture variation in adherence, completion, delays, or errors in administration, PEP completion is overall high in Madagascar and in other settings where PEP is provided at no direct costs to patients, so our estimates at 3 doses per patient vs. 4 are likely encompassing much of this variation at the annual level.
We find that estimates of vial demand based on uncorrected bite patient numbers are generally lower than the true number of vials provisioned for those clinics with substantial underreporting (Figure S2.1B). Correcting patient numbers by under-submission of forms based on the 15 day consecutive day cut-off results in estimates of vial demand closer to the observed vials provisioned.
Figure S1.1 correcting
Figure S1.2 vial ests
Excluding Category I exposures
We also excluded Category I exposures because they do not contribute to rabies risk and more importantly follow a different travel time distribution. Because they are usually actively referred to clinics (Moramanga data fig). They also present clustered on the same dates (Moramanga data fig). So we excluded by this clustering.
Figure S1.3 contacts
Figure S1.4 ttime differences
Because there was variation in estimates of reporting, and in our resulting incidence estimates we looked at how this impacted estimates of travel time effects (Supplementary Materials S4).
S3. Candidate models of bite incidence
We compared alternative model structures for predicting reported bites, specifically by including human population size as a predictor of reported bites (addPop = population size as additional covariate, onlyPop = population as only covariate, flatPop = population as offset in model), as well as by intercept type ( random = random intercept by catchment, fixed = a single fixed intercept was estimated) for the national data (see main Methods for specific formulations). We used priors centered around zero (norm(0, 10-6)) except for the precision parameter for which we used an uniform distribution (unif(0, 100)).
The Moramanga data were not fit to models with a random catchment effect, as these were data from a single catchment. We calculated Deviance Information Criterion (DIC, a metric of model fit to data) for each candidate model as well as potential scale reduction factors (psrf, a metric of model convergence) for each covariate in the model and the multivariate psrf for the whole model (cite Gelman book). All models and posterior estimates of parameters converged (Fig S3.1, both psrf and mpsrf approach 1). Population size was the poorest predictor of bites (Fig S3.2, Table S3.1) and adding population size to the model did not significantly improve the fit to the national level data (Figure S3.2), although these models did have lower DIC values for the Moramanga data (Table S3.1). However, when used to predict data that were not used to fit the models (i.e. model estimates from fitting to the Moramanga data were used to predict the national data and vice versa), only models with incidence offset by population size (flatPop) rather than scaling with population size (addPop, onlyPop) generated reasonable predictions to the data (Figure S3.3).
| Dataset | Scale | Pop effect | Intercept type | DIC |
|---|---|---|---|---|
| Moramanga | Commune | addPop | fixed | 6.42 |
| Moramanga | Commune | flatPop | fixed | 8.56 |
| Moramanga | Commune | onlyPop | fixed | 9.54 |
| National | Commune | flatPop | random | 40.75 |
| National | Commune | addPop | random | 48.62 |
| National | District | flatPop | random | 56.98 |
| National | District | addPop | random | 58.89 |
| National | Commune | flatPop | fixed | 86.17 |
| National | District | flatPop | fixed | 110.32 |
| National | Commune | addPop | fixed | 114.98 |
| National | District | onlyPop | random | 119.58 |
| National | District | addPop | fixed | 120.34 |
| National | Commune | onlyPop | random | 137.75 |
| National | District | onlyPop | fixed | 181.61 |
| National | Commune | onlyPop | fixed | 205.70 |
Figure S3.1. Potential scale reduction factors (psrf) for all models. The x-axis is the type of psrf, either multivariate which assesses the convergence of the model, or the individual which assesses the convergence for each parameter estimated from the model. Columns are by the type of model intercept (either a fixed intercept or a random intercept by catchment) and rows are the type of model structure with respect to the population covariate (addPop = population size as additional covariate, onlyPop = population as only covariate, flatPop = population as offset in model). Colors show which data set was used for fitting and the scale of the model (Moramanga = Moramanga data, Commune = National data with covariates at the commune level, District = National data with covariates at the district level).
Figure S3.2. Prediction to data used to fit each model. Log of the observed bites against the predicted bites using the mean of the posterior distributions for each parameter. Columns are by the type of model intercept (either a fixed intercept or a random intercept by catchment) and rows are the type of model structure with respect to the population covariate (addPop = population size as additional covariate, onlyPop = population as only covariate, flatPop = population as offset in model). Colors show which data set was used for fitting and the scale of the model (Moramanga = Moramanga data, Commune = National data with covariates at the commune level, District = National data with covariates at the district level).
Figure S3.3. Out of fit predictions to data. Log of the observed bites against the log of predicted bites generated from the mean of the posterior distributions from each model for data not used to fit the model. A) Estimates from commune and district model fitted to the national data predicting the Moramanga data. B) Estimates from models fitted to the Moramanga data predicting the national data at the commune and district scale. Columns are by the type of model intercept (either a fixed intercept or a random intercept by catchment) and rows are the type of model structure with respect to the population covariate (addPop = population size as additional covariate, onlyPop = population as only covariate, flatPop = population as offset in model). Colors show which data set was used for fitting (Commune = National data with covariates at the commune level, District = National data with covariates at the district level, A) and the scale of the model at which covariates were applied (Commune = covariates at the commune scale, District = covariates at the District scale).
S4. Sensitivity of these models to reporting/contact cut_offs
For all resulting data sets used to fit models, the commune scale model with population as an offset and a random catchment effect had the lowest DIC.
S5. Incrementally adding clinics (where they were added, shifts in ttimes + catchment pops, shifts in burden/reporting etc. –answer q here: why are there still abt 100 deaths)
S6. Sensitivity of burden estimates to parameter assumptions (baseline + incremental)
REFERENCES
1. Weiss DJ, Nelson A, Gibson HS, Temperley W, Peedell S, Lieber A, et al. A global map of travel time to cities to assess inequalities in accessibility in 2015. Nature. 2018;553: 333–336. doi:10.1038/nature25181
2. Linard C, Gilbert M, Snow RW, Noor AM, Tatem AJ. Population Distribution, Settlement Patterns and Accessibility across Africa in 2010. PLoS ONE. 2012.
3. Rajeev M, Edosoa G, Hanitriniaina C, Andriamandimby SF, Guis H, Ramiandrasoa R, et al. Healthcare utilization, provisioning of post-exposure prophylaxis, and estimation of human rabies burden in Madagascar. Vaccine. 2018.
4. Consortium WRM. The potential impact of improved provision of rabies post-exposure prophylaxis in Gavi-eligible countries: A modelling study. The Lancet Infectious diseases. 2019.
5. Knobel DL, Cleaveland S, Coleman PG, F‘evre EM, Meltzer MI, Miranda MEG, et al. Re-evaluating the burden of rabies in Africa and Asia. Bulletin of the World Health Organization. 2005.
6. Sambo M, Hampson K, Changalucha J, Cleaveland S, Lembo T, Lushasi K, et al. Estimating the size of dog populations in Tanzania to inform rabies control. Veterinary sciences. 2018;5: 77.
7. Changalucha J, Steenson R, Grieve E, Cleaveland S, Lembo T, Lushasi K, et al. The need to improve access to rabies post-exposure vaccines: Lessons from Tanzania. Vaccine. 2018.
8. Tarantola A, Ly S, Chan M, In S, Peng Y, Hing C, et al. Intradermal rabies post-exposure prophylaxis can be abridged with no measurable impact on clinical outcome in Cambodia, 2003-2014. Vaccine. 2018.
9. Abbas SS, Venkataramanan V, Pathak G, Kakkar M. Rabies control initiative in Tamil Nadu, India: A test case for the “One Health” approach. International Health. 2011.